/**
 * @file rad_rcslw.cc
 * \brief Source file for child class rad_rcslw
 */

#include <algorithm>
#include <cmath>
#include <fstream>
#include <iostream>
#include <cstdlib>    // exit
#include "multilinear_interpolation.h"
#include "rad_rcslw.h"

using namespace std;

#define GETSTRINGEDVALUE(themacro) STRINGIFY(themacro)
#define STRINGIFY(thestring) #thestring

////////////////////////////////////////////////////////////////////////////////
///
/// Constructor function
/// @param p_nGG        \input   number of gray gases not including the gray gas
/// @param TbTref       \input   black temperature and reference temperature (K); set once, assumed constant
/// @param p_P          \input   system pressure (Pa); set once, assumed constant
/// @param fvsoot       \input   soot volume fraction
/// @param xH2O         \input   H2O mole fraction
/// @param xCO2         \input   CO2 mole fraction
/// @param xCO          \input   CO mole fraction
///
////////////////////////////////////////////////////////////////////////////////

rad_rcslw::rad_rcslw(const int    p_nGG,
                     const double TbTref,
                     const double p_P,
                     const double fvsoot,
                     const double xH2O,
                     const double xCO2,
                     const double xCO) :
        rad(p_nGG, p_nGG + 1){

    P     = p_P/101325;      // atm
    Tref  = TbTref;
    Tb    = TbTref;

    Cmin = 0.0001;
    Cmax = 1000.0;

    P_table  = vector<double>{0.1, 0.25, 0.5, 1,2,4,8,15,30,50};
    C_table  = vector<double>(71);
    for(int  i=0; i<71; i++)
        C_table[i] = Cmin*pow(Cmax/Cmin, i/70.0);

    Tg_table = vector<double>(28);
    for(int i=0; i<28; i++)
        Tg_table[i] = 300 + i*(3000.0-300.0)/27;

    Tb_table = vector<double>(28);
    for(int i=0; i<28; i++)
        Tb_table[i] = 300 + i*(3000.0-300.0)/27;

    xH2O_table = vector<double>{0, 0.05, 0.1, 0.2, 0.3, 0.4, 0.6, 0.8, 1.0};

    nP     =  P_table.size();        // 10 
    nC     =  C_table.size();        // 71
    nTg    =  Tg_table.size();       // 28
    nTb    =  Tb_table.size();       // 28
    ny_H2O =  xH2O_table.size();     //  9

    set_Falbdf_CO2_CO_H2O_at_P();

    Fmin = get_F_albdf(Cmin, Tref, Tref, xCO2, xCO, xH2O, fvsoot);
    Fmax = get_F_albdf(Cmax, Tref, Tref, xCO2, xCO, xH2O, fvsoot);

    set_Fpts();

}

////////////////////////////////////////////////////////////////////////////////
///
/// **This is the class interface function**
/// Given the gas state, set the k and a vectors.
///
/// These can then be accessed by the user.
///
/// Return through arg list the local gray gas coefficients (kabs) and the local weights (awts).
/// @param kabs            \output absorption coefficient (1/m) for band/gas iband: ranges from 0 to nGG inclusive
/// @param awts            \output weight (unitless; total sums to 1) for band/gas iband: ranges from 0 to nGG inclusive
/// @param iband           \input which band to compute
/// @param T               \input  gas temperature
/// @param P_not_used      \input  Pressure (Pa)      NOT USED; HERE FOR INTERFACE; P IS SET BY CONSTRUCTOR
/// @param fvsoot          \input  soot volume fraction = rho*Ysoot/rhosoot
/// @param xH2O            \input  mole fraction H2O
/// @param xCO2            \input  mole fraction CO2
/// @param xCO             \input  mole fraction CO
/// @param xCH4_not_used   \input  mole fraction CH4  NOT USED; HERE FOR INTERFACE; (... pass in 0.0)
///
////////////////////////////////////////////////////////////////////////////////

void rad_rcslw::get_k_a_oneband(double         &kabs,
                                double         &awts,
                                const int      iband,
                                const double   T,
                                const double   P_not_used,
                                const double   fvsoot,
                                const double   xH2O,
                                const double   xCO2,
                                const double   xCO,
                                const double   xCH4_not_used){

    if(iband < 0 || iband >= nGGa) {
        cerr << "\n\n***** rad_rcslw::get_k_a_oneband: iband out of range *****\n" << endl; 
        exit(0); 
    }

    //--------------- kabs

    if(iband==0)
        kabs = 0.0;
    else{
        double Nconc = P*101325/8.31446/T;    // mol/m3
        double C = get_FI_albdf(F_pts[iband-1], T, Tb, xCO2, xCO, xH2O, fvsoot);
        kabs = Nconc * C;
    }

    //--------------- awts

    if(iband==0){
        double Ct = get_FI_albdf(Ft_pts[0], T, Tb, xCO2, xCO, xH2O, fvsoot);
        awts = get_F_albdf(Ct, T, T, xCO2, xCO, xH2O, fvsoot);
    }
    else{
        double Ct_i   = get_FI_albdf(Ft_pts[iband],   T, Tb, xCO2, xCO, xH2O, fvsoot);
        double Ct_im1 = get_FI_albdf(Ft_pts[iband-1], T, Tb, xCO2, xCO, xH2O, fvsoot);
        awts = get_F_albdf(Ct_i,   T, T, xCO2, xCO, xH2O, fvsoot) - 
               get_F_albdf(Ct_im1, T, T, xCO2, xCO, xH2O, fvsoot);
    }

}

////////////////////////////////////////////////////////////////////////////////
///
/// **This is the class interface function**
/// Given the gas state, set the k and a vectors.
///
/// These can then be accessed by the user.
///
/// Return through arg list the local gray gas coefficients (kabs) and the local weights (awts).
/// @param kabs            \output absorption coefficients (1/m) for nGG+1 (nGG gray gases + clear gas)
/// @param awts            \output weights (unitless; sums to 1) for nGG+1 (nGG gray gases + clear gas)
/// @param T               \input  gas temperature
/// @param P_not_used      \input  Pressure (Pa)      NOT USED; HERE FOR INTERFACE; P IS SET BY CONSTRUCTOR
/// @param fvsoot          \input  soot volume fraction = rho*Ysoot/rhosoot
/// @param xH2O            \input  mole fraction H2O
/// @param xCO2            \input  mole fraction CO2
/// @param xCO             \input  mole fraction CO
/// @param xCH4_not_used   \input  mole fraction CH4  NOT USED; HERE FOR INTERFACE; (... pass in 0.0)
///
////////////////////////////////////////////////////////////////////////////////

void rad_rcslw::get_k_a(vector<double> &kabs,
                        vector<double> &awts,
                        const double   T,
                        const double   P_not_used,
                        const double   fvsoot,
                        const double   xH2O,
                        const double   xCO2,
                        const double   xCO,
                        const double   xCH4_not_used){

    vector<double> C(nGG);
    vector<double> Ct(nGGa);
    vector<double> FCt(nGGa);

    for(int j=0; j < nGG; j++)
        C[j] = get_FI_albdf(F_pts[j], T, Tb, xCO2, xCO, xH2O, fvsoot);
    
    for(int j=0; j < nGGa; j++)
        Ct[j] = get_FI_albdf(Ft_pts[j], T, Tb, xCO2, xCO, xH2O, fvsoot);
    
    double Nconc = P*101325/8.31446/T;    // mol/m3
    
    kabs.resize(nGGa);
    kabs[0] = 0.0;
    for(int j=1; j < nGGa; j++)
        kabs[j] = Nconc * C[j-1];

    for(int j=0; j < nGGa; j++)
        FCt[j] = get_F_albdf(Ct[j], T, T, xCO2, xCO, xH2O, fvsoot);

    awts.resize(nGGa);
    awts[0] = FCt[0];
    for(int j=1; j < nGGa; j++)
        awts[j] = FCt[j] - FCt[j-1];

    //------------------ contents below could be used instead for this whole function, but slower

    //kabs.resize(nGGa);
    //awts.resize(nGGa);

    //double k, a;
    //for(int i=0; i<nGGa; i++){
    //    get_k_a_oneband(k, a, i, T, P, fvsoot, xH2O, xCO2, xCO, xCH4_not_used);
    //    kabs[i] = k;
    //    awts[i] = a;
    //}
}

/////////////////////////////////////////////////////////////
///
/// Get the albdf function F
/// @param C      \input  cross section (m2/mol)
/// @param Tg     \input  gas temperature  
/// @param Tb     \input  black temperature
/// @param xCO2   \input  mole fraction CO2
/// @param xCO    \input  mole fraction CO
/// @param xH2O   \input  mole fraction H2O
/// @param fvsoot \input  soot volume fraction
/// @return the albdf function F
///
/////////////////////////////////////////////////////////////

double rad_rcslw::get_F_albdf(const double C, double Tg, double Tb,
                              double xCO2, double xCO, double xH2O,
                              double fvsoot){


    if (xCO2 <= 1E-20)               xCO2 = 1E-20;
    if (xCO <= 1e-20)                xCO = 1e-20;
    if (xH2O <= 1e-20)               xH2O = 1e-20;
    if (fvsoot < 0.0)                fvsoot = 0.0;

    if (Tg < Tg_table[0])            Tg = Tg_table[0];
    if (Tg > Tg_table.back())        Tg = Tg_table.back();

    if (Tb < Tb_table[0])            Tb = Tb_table[0];
    if (Tb > Tb_table.back())        Tb = Tb_table.back();

    if (xH2O < xH2O_table[0])        xH2O = xH2O_table[0];
    if (xH2O > xH2O_table.back()) xH2O = xH2O_table.back();

    double CxCO2 = C/xCO2;
    double CxCO  = C/xCO;
    double CxH2O = C/xH2O;

    if (CxCO2  < C_table[0])       CxCO2 = C_table[0];
    if (CxCO2  > C_table.back())   CxCO2 = C_table.back();
    if (CxCO   < C_table[0])       CxCO  = C_table[0];
    if (CxCO   > C_table.back())   CxCO  = C_table.back();
    if (CxH2O  < C_table[0])       CxH2O = C_table[0];
    if (CxH2O  > C_table.back())   CxH2O = C_table.back();

    double F_CO2, F_CO, F_H2O;

    F_CO2 = LI_3D(nTg, nTb, nC, &Tg_table[0], &Tb_table[0], &C_table[0], &Falbdf_CO2[0],        Tg, Tb, CxCO2);
    F_CO  = LI_3D(nTg, nTb, nC, &Tg_table[0], &Tb_table[0], &C_table[0], &Falbdf_CO[0],         Tg, Tb, CxCO);
    F_H2O = LI_4D(ny_H2O, nTg, nTb, nC, &xH2O_table[0], &Tg_table[0], &Tb_table[0], &C_table[0], &Falbdf_H2O[0],  xH2O, Tg, Tb, CxH2O);

    return F_CO2 * F_CO * F_H2O * F_albdf_soot(C, Tg, Tb, fvsoot);
}


///////////////////////////////////////////////////////
///
/// Inverse F_albdf: pass in albdf F and get out Cross section C
/// @param F       \input   albdf
/// @param Tg      \input   gas temperature
/// @param Tb      \input   black temperature
/// @param xCO2    \input   mole fraction CO2
/// @param xCO     \input   mole fraction CO
/// @param xH2O    \input   mole fraction H2O
/// @param fvsoot  \input   soot volume fraction = rho*Ysoot/rhoSoot
/// @return Cross section C (m2/mol).
///
///////////////////////////////////////////////////////

double rad_rcslw::get_FI_albdf(const double F, const double Tg, const double Tb,
                               const double xCO2, const double xCO, const double xH2O,
                               const double fvsoot){

   if (F <= get_F_albdf(Cmin, Tg, Tb, xCO2, xCO, xH2O, fvsoot))
       return(Cmin);
   if (F >= get_F_albdf(Cmax, Tg, Tb, xCO2, xCO, xH2O, fvsoot))
        return(Cmax);

   // Find table location

   int iLo = 0;
   int iHi = nC - 1;      

   int maxit = 100;
   int it;
   for (it=0; it < maxit; it++){
       if (iHi-iLo == 1)
           break;
       int iMi = (iLo+iHi)/2;
       if (F > get_F_albdf(C_table[iMi], Tg, Tb, xCO2, xCO, xH2O, fvsoot))
           iLo = iMi;
       else
           iHi = iMi;
   }

   int iMi = (iLo+iHi)/2;
   if(it == maxit){
       cout << "WARNING, NO CONVERGENCE IN" << maxit << "iterations" << "\n";
       return C_table[iMi] ;
   }

   // Now interpolate to the solution

   double Flo = get_F_albdf(C_table[iLo], Tg, Tb, xCO2, xCO, xH2O, fvsoot);
   double Fhi = get_F_albdf(C_table[iHi], Tg, Tb, xCO2, xCO, xH2O, fvsoot);
   double cc = C_table[iLo] + (F-Flo) * (C_table[iHi]-C_table[iLo])/(Fhi-Flo);

   // cout << "relative error:" << (F - get_F_albdf(cc, Tg, Tb, xCO2, xCO, xH2O, fvsoot))/F << "\n";
   // return cc; //doldb

   //----------- but F(cc) is not equal to F! So give it another interpolation:
   // rel error ~ 0.1% --> 1E-6

   double FF = get_F_albdf(cc, Tg, Tb, xCO2, xCO, xH2O, fvsoot);
   double ccc = C_table[iLo] + (F-Flo)*(cc-C_table[iLo])/(FF-Flo);

   return ccc;
   
   ////----------- and one more for fun...
   //// rel error ~ 1E-6 --> 1E-10

   //double FFF = get_F_albdf(ccc, Tg, Tb, xCO2, xCO, xH2O, fvsoot);
   //double cccc = cc + (F-FF)*(ccc-cc)/(FFF-FF);

   ////cout << "relative error:" << (F - get_F_albdf(ccc, Tg, Tb, xCO2, xCO, xH2O, fvsoot))/F << "\n";

   //return cccc;

}

////////////////////////////////////////////////////////////////
///
/// Set grid of F points based on Gauss-Legendre Quadrature.
///
////////////////////////////////////////////////////////////////

void rad_rcslw::set_Fpts(){

    if(nGG > 25){
        cout << endl << "ERROR: nGG too high for quadrature: nGG max is currently 25" << endl;
        exit(0);
    }

    vector<vector<double> > W(26);
    
    W[0]  = vector<double> {0.0,0.0};
    W[1]  = vector<double> {1.0000000000000000E+00, 1.0000000000000000E+00};
    W[2]  = vector<double> {3.4785484513745368E-01, 6.5214515486254621E-01, 6.5214515486254621E-01, 3.4785484513745368E-01};
    W[3]  = vector<double> {1.7132449237916975E-01, 3.6076157304813894E-01, 4.6791393457269137E-01, 4.6791393457269137E-01, 3.6076157304813894E-01, 1.7132449237916975E-01};
    W[4]  = vector<double> {1.0122853629037669E-01, 2.2238103445337434E-01, 3.1370664587788705E-01, 3.6268378337836177E-01, 3.6268378337836177E-01, 3.1370664587788705E-01, 2.2238103445337434E-01, 1.0122853629037669E-01};
    W[5]  = vector<double> {6.6671344308688069E-02, 1.4945134915058036E-01, 2.1908636251598201E-01, 2.6926671930999652E-01, 2.9552422471475298E-01, 2.9552422471475298E-01, 2.6926671930999652E-01, 2.1908636251598201E-01, 1.4945134915058036E-01, 6.6671344308688069E-02}; 
    W[6]  = vector<double> {4.7175336386512022E-02, 1.0693932599531888E-01, 1.6007832854334611E-01, 2.0316742672306565E-01, 2.3349253653835464E-01, 2.4914704581340269E-01, 2.4914704581340269E-01, 2.3349253653835464E-01, 2.0316742672306565E-01, 1.6007832854334611E-01, 1.0693932599531888E-01, 4.7175336386512022E-02};
    W[7]  = vector<double> {3.5119460331752374E-02, 8.0158087159760305E-02, 1.2151857068790296E-01, 1.5720316715819341E-01, 1.8553839747793763E-01, 2.0519846372129555E-01, 2.1526385346315766E-01, 2.1526385346315766E-01, 2.0519846372129555E-01, 1.8553839747793763E-01, 1.5720316715819341E-01, 1.2151857068790296E-01, 8.0158087159760305E-02, 3.5119460331752374E-02};
    W[8]  = vector<double> {2.7152459411754037E-02, 6.2253523938647706E-02, 9.5158511682492591E-02, 1.2462897125553403E-01, 1.4959598881657676E-01, 1.6915651939500262E-01, 1.8260341504492361E-01, 1.8945061045506859E-01, 1.8945061045506859E-01, 1.8260341504492361E-01, 1.6915651939500262E-01, 1.4959598881657676E-01, 1.2462897125553403E-01, 9.5158511682492591E-02, 6.2253523938647706E-02, 2.7152459411754037E-02};
    W[9]  = vector<double> {2.1616013526484130E-02, 4.9714548894969221E-02, 7.6425730254889246E-02, 1.0094204410628699E-01, 1.2255520671147836E-01, 1.4064291467065063E-01, 1.5468467512626521E-01, 1.6427648374583273E-01, 1.6914238296314363E-01, 1.6914238296314363E-01, 1.6427648374583273E-01, 1.5468467512626521E-01, 1.4064291467065063E-01, 1.2255520671147836E-01, 1.0094204410628699E-01, 7.6425730254889246E-02, 4.9714548894969221E-02, 2.1616013526484130E-02};
    W[10] = vector<double> {1.7614007139153273E-02, 4.0601429800386217E-02, 6.2672048334109443E-02, 8.3276741576704671E-02, 1.0193011981724026E-01, 1.1819453196151825E-01, 1.3168863844917653E-01, 1.4209610931838187E-01, 1.4917298647260366E-01, 1.5275338713072578E-01, 1.5275338713072578E-01, 1.4917298647260366E-01, 1.4209610931838187E-01, 1.3168863844917653E-01, 1.1819453196151825E-01, 1.0193011981724026E-01, 8.3276741576704671E-02, 6.2672048334109443E-02, 4.0601429800386217E-02, 1.7614007139153273E-02}; 
    W[11] = vector<double> {1.4627995298274705E-02, 3.3774901584815178E-02, 5.2293335152682870E-02, 6.9796468424520197E-02, 8.5941606217067396E-02, 1.0041414444288072E-01, 1.1293229608053883E-01, 1.2325237681051199E-01, 1.3117350478706188E-01, 1.3654149834601478E-01, 1.3925187285563156E-01, 1.3925187285563156E-01, 1.3654149834601478E-01, 1.3117350478706188E-01, 1.2325237681051199E-01, 1.1293229608053883E-01, 1.0041414444288072E-01, 8.5941606217067396E-02, 6.9796468424520197E-02, 5.2293335152682870E-02, 3.3774901584815178E-02, 1.4627995298274705E-02}; 
    W[12] = vector<double> {1.2341229799987091E-02, 2.8531388628933743E-02, 4.4277438817419551E-02, 5.9298584915436742E-02, 7.3346481411080411E-02, 8.6190161531953288E-02, 9.7618652104114065E-02, 1.0744427011596561E-01, 1.1550566805372561E-01, 1.2167047292780342E-01, 1.2583745634682830E-01, 1.2793819534675221E-01, 1.2793819534675221E-01, 1.2583745634682830E-01, 1.2167047292780342E-01, 1.1550566805372561E-01, 1.0744427011596561E-01, 9.7618652104114065E-02, 8.6190161531953288E-02, 7.3346481411080411E-02, 5.9298584915436742E-02, 4.4277438817419551E-02, 2.8531388628933743E-02, 1.2341229799987091E-02}; 
    W[13] = vector<double> {1.0551372617343395E-02, 2.4417851092631938E-02, 3.7962383294363120E-02, 5.0975825297148079E-02, 6.3274046329574674E-02, 7.4684149765659763E-02, 8.5045894313485068E-02, 9.4213800355914160E-02, 1.0205916109442532E-01, 1.0847184052857647E-01, 1.1336181654631956E-01, 1.1666044348529646E-01, 1.1832141527926213E-01, 1.1832141527926213E-01, 1.1666044348529646E-01, 1.1336181654631956E-01, 1.0847184052857647E-01, 1.0205916109442532E-01, 9.4213800355914160E-02, 8.5045894313485068E-02, 7.4684149765659763E-02, 6.3274046329574674E-02, 5.0975825297148079E-02, 3.7962383294363120E-02, 2.4417851092631938E-02, 1.0551372617343395E-02}; 
    W[14] = vector<double> {9.1242825930943974E-03, 2.1132112592771271E-02, 3.2901427782304517E-02, 4.4272934759003985E-02, 5.5107345675716936E-02, 6.5272923966999755E-02, 7.4646214234568811E-02, 8.3113417228900935E-02, 9.0571744393032852E-02, 9.6930657997929923E-02, 1.0211296757806078E-01, 1.0605576592284637E-01, 1.0871119225829413E-01, 1.1004701301647524E-01, 1.1004701301647524E-01, 1.0871119225829413E-01, 1.0605576592284637E-01, 1.0211296757806078E-01, 9.6930657997929923E-02, 9.0571744393032852E-02, 8.3113417228900935E-02, 7.4646214234568811E-02, 6.5272923966999755E-02, 5.5107345675716936E-02, 4.4272934759003985E-02, 3.2901427782304517E-02, 2.1132112592771271E-02, 9.1242825930943974E-03};
    W[15] = vector<double> {7.9681924961695228E-03, 1.8466468311091087E-02, 2.8784707883322873E-02, 3.8799192569626793E-02, 4.8402672830594434E-02, 5.7493156217619093E-02, 6.5974229882180324E-02, 7.3755974737704802E-02, 8.0755895229419811E-02, 8.6899787201082698E-02, 9.2122522237785789E-02, 9.6368737174643990E-02, 9.9593420586794934E-02, 1.0176238974840521E-01, 1.0285265289355848E-01, 1.0285265289355848E-01, 1.0176238974840521E-01, 9.9593420586794934E-02, 9.6368737174643990E-02, 9.2122522237785789E-02, 8.6899787201082698E-02, 8.0755895229419811E-02, 7.3755974737704802E-02, 6.5974229882180324E-02, 5.7493156217619093E-02, 4.8402672830594434E-02, 3.8799192569626793E-02, 2.8784707883322873E-02, 1.8466468311091087E-02, 7.9681924961695228E-03}; 
    W[16] = vector<double> {7.0186100094692984E-03, 1.6274394730905965E-02, 2.5392065309262427E-02, 3.4273862913021626E-02, 4.2835898022226426E-02, 5.0998059262376244E-02, 5.8684093478535704E-02, 6.5822222776361752E-02, 7.2345794108848449E-02, 7.8193895787070311E-02, 8.3311924226946846E-02, 8.7652093004403908E-02, 9.1173878695763863E-02, 9.3844399080804566E-02, 9.5638720079274833E-02, 9.6540088514727812E-02, 9.6540088514727812E-02, 9.5638720079274833E-02, 9.3844399080804566E-02, 9.1173878695763863E-02, 8.7652093004403908E-02, 8.3311924226946846E-02, 7.8193895787070311E-02, 7.2345794108848449E-02, 6.5822222776361752E-02, 5.8684093478535704E-02, 5.0998059262376244E-02, 4.2835898022226426E-02, 3.4273862913021626E-02, 2.5392065309262427E-02, 1.6274394730905965E-02, 7.0186100094692984E-03};
    W[17] = vector<double> {6.2291405559100326E-03, 1.4450162748594548E-02, 2.2563721985495038E-02, 3.0491380638445777E-02, 3.8166593796386906E-02, 4.5525611523353778E-02, 5.2507414572678060E-02, 5.9054135827524654E-02, 6.5111521554076457E-02, 7.0629375814255727E-02, 7.5561974660031797E-02, 7.9868444339771819E-02, 8.3513099699845522E-02, 8.6465739747035669E-02, 8.8701897835693738E-02, 9.0203044370640612E-02, 9.0956740330259786E-02, 9.0956740330259786E-02, 9.0203044370640612E-02, 8.8701897835693738E-02, 8.6465739747035669E-02, 8.3513099699845522E-02, 7.9868444339771819E-02, 7.5561974660031797E-02, 7.0629375814255727E-02, 6.5111521554076457E-02, 5.9054135827524654E-02, 5.2507414572678060E-02, 4.5525611523353778E-02, 3.8166593796386906E-02, 3.0491380638445777E-02, 2.2563721985495038E-02, 1.4450162748594548E-02, 6.2291405559100326E-03}; 
    W[18] = vector<double> {5.5657196642477837E-03, 1.2915947284064104E-02, 2.0181515297735174E-02, 2.7298621498568355E-02, 3.4213810770307482E-02, 4.0875750923645232E-02, 4.7235083490266047E-02, 5.3244713977759678E-02, 5.8860144245324549E-02, 6.4039797355015429E-02, 6.8745323835736297E-02, 7.2941885005653004E-02, 7.6598410645870627E-02, 7.9687828912071559E-02, 8.2187266704339651E-02, 8.4078218979661792E-02, 8.5346685739338499E-02, 8.5983275670394627E-02, 8.5983275670394627E-02, 8.5346685739338499E-02, 8.4078218979661792E-02, 8.2187266704339651E-02, 7.9687828912071559E-02, 7.6598410645870627E-02, 7.2941885005653004E-02, 6.8745323835736297E-02, 6.4039797355015429E-02, 5.8860144245324549E-02, 5.3244713977759678E-02, 4.7235083490266047E-02, 4.0875750923645232E-02, 3.4213810770307482E-02, 2.7298621498568355E-02, 2.0181515297735174E-02, 1.2915947284064104E-02, 5.5657196642477837E-03}; 
    W[19] = vector<double> {5.0028807496389329E-03, 1.1613444716468455E-02, 1.8156577709613410E-02, 2.4579739738232007E-02, 3.0839500545175719E-02, 3.6894081594024859E-02, 4.2703158504674703E-02, 4.8228061860758530E-02, 5.3432019910332196E-02, 5.8280399146997154E-02, 6.2740933392133172E-02, 6.6783937979140381E-02, 7.0382507066898942E-02, 7.3512692584743411E-02, 7.6153663548446479E-02, 7.8287844658210981E-02, 7.9901033243527819E-02, 8.0982493770597061E-02, 8.1525029280385797E-02, 8.1525029280385797E-02, 8.0982493770597061E-02, 7.9901033243527819E-02, 7.8287844658210981E-02, 7.6153663548446479E-02, 7.3512692584743411E-02, 7.0382507066898942E-02, 6.6783937979140381E-02, 6.2740933392133172E-02, 5.8280399146997154E-02, 5.3432019910332196E-02, 4.8228061860758530E-02, 4.2703158504674703E-02, 3.6894081594024859E-02, 3.0839500545175719E-02, 2.4579739738232007E-02, 1.8156577709613410E-02, 1.1613444716468455E-02, 5.0028807496389329E-03};
    W[20] = vector<double> {4.5212770985300181E-03, 1.0498284531151609E-02, 1.6421058381907345E-02, 2.2245849194166653E-02, 2.7937006980023528E-02, 3.3460195282547678E-02, 3.8782167974472377E-02, 4.3870908185673324E-02, 4.8695807635072405E-02, 5.3227846983937115E-02, 5.7439769099391892E-02, 6.1306242492929319E-02, 6.4804013456601486E-02, 6.7912045815234398E-02, 7.0611647391287169E-02, 7.2886582395804478E-02, 7.4723169057968677E-02, 7.6110361900626741E-02, 7.7039818164248389E-02, 7.7505947978425332E-02, 7.7505947978425332E-02, 7.7039818164248389E-02, 7.6110361900626741E-02, 7.4723169057968677E-02, 7.2886582395804478E-02, 7.0611647391287169E-02, 6.7912045815234398E-02, 6.4804013456601486E-02, 6.1306242492929319E-02, 5.7439769099391892E-02, 5.3227846983937115E-02, 4.8695807635072405E-02, 4.3870908185673324E-02, 3.8782167974472377E-02, 3.3460195282547678E-02, 2.7937006980023528E-02, 2.2245849194166653E-02, 1.6421058381907345E-02, 1.0498284531151609E-02, 4.5212770985300181E-03};
    
    W[21] = vector<double> {4.105998604646913e-03, 9.536220301748407e-03, 1.4922443697357294e-02, 2.0227869569052214e-02, 2.5422959526113627e-02, 3.0479240699603335e-02, 3.5369071097592e-02, 4.0065735180692286e-02, 4.454357777196598e-02, 4.877814079280347e-02, 5.274629569917415e-02, 5.642636935801852e-02, 5.979826222758681e-02, 6.284355804500275e-02, 6.554562436490913e-02, 6.788970337652217e-02, 6.986299249259438e-02, 7.14547142651712e-02, 7.265617524380438e-02, 7.346081345346782e-02, 7.386423423217317e-02, 7.386423423217317e-02, 7.346081345346782e-02, 7.265617524380438e-02, 7.14547142651712e-02, 6.986299249259438e-02, 6.788970337652217e-02, 6.554562436490913e-02, 6.284355804500275e-02, 5.979826222758681e-02, 5.642636935801852e-02, 5.274629569917415e-02, 4.877814079280347e-02, 4.454357777196598e-02, 4.0065735180692286e-02, 3.5369071097592e-02, 3.0479240699603335e-02, 2.5422959526113627e-02, 2.0227869569052214e-02, 1.4922443697357294e-02, 9.536220301748407e-03, 4.105998604646913e-03};
    W[22] = vector<double> {3.7454048031164908e-03, 8.700481367523676e-03, 1.3619586755579383e-02, 1.8471481736814985e-02, 2.323148190201923e-02, 2.787578282128097e-02, 3.238122281206993e-02, 3.672534781380876e-02, 4.0886512310346054e-02, 4.484398408197005e-02, 4.857804644835181e-02, 5.207009609170433e-02, 5.53027355637279e-02, 5.825985987759541e-02, 6.092673670156192e-02, 6.329007973320361e-02, 6.533811487918127e-02, 6.706063890629348e-02, 6.844907026936646e-02, 6.949649186157239e-02, 7.019768547355804e-02, 7.054915778935386e-02, 7.054915778935386e-02, 7.019768547355804e-02, 6.949649186157239e-02, 6.844907026936646e-02, 6.706063890629348e-02, 6.533811487918127e-02, 6.329007973320361e-02, 6.092673670156192e-02, 5.825985987759541e-02, 5.53027355637279e-02, 5.207009609170433e-02, 4.857804644835181e-02, 4.484398408197005e-02, 4.0886512310346054e-02, 3.672534781380876e-02, 3.238122281206993e-02, 2.787578282128097e-02, 2.323148190201923e-02, 1.8471481736814985e-02, 1.3619586755579383e-02, 8.700481367523676e-03, 3.7454048031164908e-03};
    W[23] = vector<double> {3.4303008681091456e-03, 7.969898229724492e-03, 1.2479883770988784e-02, 1.6933514007836395e-02, 2.1309998754136483e-02, 2.5589286397129797e-02, 2.975182955220257e-02, 3.37786279991068e-02, 3.765130535738592e-02, 4.135219010967875e-02, 4.486439527731805e-02, 4.817189510171219e-02, 5.1259598007142824e-02, 5.411341538585653e-02, 5.6720325843991094e-02, 5.906843459554615e-02, 6.1147027724650325e-02, 6.294662106439443e-02, 6.445900346713893e-02, 6.567727426778112e-02, 6.659587476845476e-02, 6.721061360067808e-02, 6.751868584903632e-02, 6.751868584903632e-02, 6.721061360067808e-02, 6.659587476845476e-02, 6.567727426778112e-02, 6.445900346713893e-02, 6.294662106439443e-02, 6.1147027724650325e-02, 5.906843459554615e-02, 5.6720325843991094e-02, 5.411341538585653e-02, 5.1259598007142824e-02, 4.817189510171219e-02, 4.486439527731805e-02, 4.135219010967875e-02, 3.765130535738592e-02, 3.37786279991068e-02, 2.975182955220257e-02, 2.5589286397129797e-02, 2.1309998754136483e-02, 1.6933514007836395e-02, 1.2479883770988784e-02, 7.969898229724492e-03, 3.4303008681091456e-03};
    W[24] = vector<double> {3.1533460523091796e-03, 7.327553901276492e-03, 1.1477234579234974e-02, 1.5579315722942928e-02, 1.9616160457355297e-02, 2.3570760839324092e-02, 2.7426509708356882e-02, 3.116722783279834e-02, 3.477722256477066e-02, 3.8241351065830674e-02, 4.1545082943464554e-02, 4.46745608566941e-02, 4.7616658492490284e-02, 5.035903555385428e-02, 5.289018948519349e-02, 5.5199503699984054e-02, 5.727729210040293e-02, 5.911483969839548e-02, 6.070443916589358e-02, 6.2039423159892464e-02, 6.311419228625378e-02, 6.392423858464795e-02, 6.446616443594984e-02, 6.473769681268368e-02, 6.473769681268368e-02, 6.446616443594984e-02, 6.392423858464795e-02, 6.311419228625378e-02, 6.2039423159892464e-02, 6.070443916589358e-02, 5.911483969839548e-02, 5.727729210040293e-02, 5.5199503699984054e-02, 5.289018948519349e-02, 5.035903555385428e-02, 4.7616658492490284e-02, 4.46745608566941e-02, 4.1545082943464554e-02, 3.8241351065830674e-02, 3.477722256477066e-02, 3.116722783279834e-02, 2.7426509708356882e-02, 2.3570760839324092e-02, 1.9616160457355297e-02, 1.5579315722942928e-02, 1.1477234579234974e-02, 7.327553901276492e-03, 3.1533460523091796e-03};
    W[25] = vector<double> {2.9086225531578685e-03, 6.759799195745505e-03, 1.059054838365167e-02, 1.438082276148571e-02, 1.811556071348957e-02, 2.1780243170124582e-02, 2.536067357001259e-02, 2.8842993580534923e-02, 3.2213728223577896e-02, 3.5459835615146026e-02, 3.8568756612587844e-02, 4.15284630901474e-02, 4.43275043388031e-02, 4.695505130394827e-02, 4.9400938449466136e-02, 5.165570306958092e-02, 5.371062188899592e-02, 5.555774480621228e-02, 5.718992564772815e-02, 5.8600849813222215e-02, 5.9785058704265225e-02, 6.073797084176991e-02, 6.145589959031641e-02, 6.193606742068298e-02, 6.2176616655347e-02, 6.2176616655347e-02, 6.193606742068298e-02, 6.145589959031641e-02, 6.073797084176991e-02, 5.9785058704265225e-02, 5.8600849813222215e-02, 5.718992564772815e-02, 5.555774480621228e-02, 5.371062188899592e-02, 5.165570306958092e-02, 4.9400938449466136e-02, 4.695505130394827e-02, 4.43275043388031e-02, 4.15284630901474e-02, 3.8568756612587844e-02, 3.5459835615146026e-02, 3.2213728223577896e-02, 2.8842993580534923e-02, 2.536067357001259e-02, 2.1780243170124582e-02, 1.811556071348957e-02, 1.438082276148571e-02, 1.059054838365167e-02, 6.759799195745505e-03, 2.9086225531578685e-03};

    vector<vector<double> > X(26);
    
    X[0]  = vector<double> {0.0, 0.0};
    X[1]  = vector<double> {-5.7735026918962573E-01,  5.7735026918962573E-01};
    X[2]  = vector<double> {-8.6113631159405257E-01, -3.3998104358485626E-01,  3.3998104358485626E-01,  8.6113631159405257E-01};
    X[3]  = vector<double> {-9.3246951420315205E-01, -6.6120938646626448E-01, -2.3861918608319693E-01,  2.3861918608319693E-01,  6.6120938646626448E-01,  9.3246951420315205E-01};
    X[4]  = vector<double> {-9.6028985649753618E-01, -7.9666647741362673E-01, -5.2553240991632899E-01, -1.8343464249564978E-01,  1.8343464249564978E-01,  5.2553240991632899E-01,  7.9666647741362673E-01,  9.6028985649753618E-01};
    X[5]  = vector<double> {-9.7390652851717174E-01, -8.6506336668898454E-01, -6.7940956829902444E-01, -4.3339539412924721E-01, -1.4887433898163122E-01,  1.4887433898163122E-01,  4.3339539412924721E-01,  6.7940956829902444E-01,  8.6506336668898454E-01,  9.7390652851717174E-01};
    X[6]  = vector<double> {-9.8156063424671924E-01, -9.0411725637047480E-01, -7.6990267419430469E-01, -5.8731795428661748E-01, -3.6783149899818018E-01, -1.2523340851146891E-01,  1.2523340851146891E-01,  3.6783149899818018E-01,  5.8731795428661748E-01,  7.6990267419430469E-01,  9.0411725637047480E-01,  9.8156063424671924E-01};
    X[7]  = vector<double> {-9.8628380869681231E-01, -9.2843488366357352E-01, -8.2720131506976502E-01, -6.8729290481168548E-01, -5.1524863635815410E-01, -3.1911236892788974E-01, -1.0805494870734367E-01,  1.0805494870734367E-01,  3.1911236892788974E-01,  5.1524863635815410E-01,  6.8729290481168548E-01,  8.2720131506976502E-01,  9.2843488366357352E-01,  9.8628380869681231E-01};
    X[8]  = vector<double> {-9.8940093499164994E-01, -9.4457502307323260E-01, -8.6563120238783176E-01, -7.5540440835500300E-01, -6.1787624440264377E-01, -4.5801677765722737E-01, -2.8160355077925892E-01, -9.5012509837637454E-02,  9.5012509837637454E-02,  2.8160355077925892E-01,  4.5801677765722737E-01,  6.1787624440264377E-01,  7.5540440835500300E-01,  8.6563120238783176E-01,  9.4457502307323260E-01,  9.8940093499164994E-01};
    X[9]  = vector<double> {-9.9156516842093090E-01, -9.5582394957139782E-01, -8.9260246649755570E-01, -8.0370495897252314E-01, -6.9168704306035322E-01, -5.5977083107394754E-01, -4.1175116146284263E-01, -2.5188622569150548E-01, -8.4775013041735292E-02,  8.4775013041735292E-02,  2.5188622569150548E-01,  4.1175116146284263E-01,  5.5977083107394754E-01,  6.9168704306035322E-01,  8.0370495897252314E-01,  8.9260246649755570E-01,  9.5582394957139782E-01,  9.9156516842093090E-01};
    X[10] = vector<double> {-9.9312859918509488E-01, -9.6397192727791381E-01, -9.1223442825132584E-01, -8.3911697182221878E-01, -7.4633190646015080E-01, -6.3605368072651502E-01, -5.1086700195082713E-01, -3.7370608871541955E-01, -2.2778585114164510E-01, -7.6526521133497338E-02,  7.6526521133497338E-02,  2.2778585114164510E-01,  3.7370608871541955E-01,  5.1086700195082713E-01,  6.3605368072651502E-01,  7.4633190646015080E-01,  8.3911697182221878E-01,  9.1223442825132584E-01,  9.6397192727791381E-01,  9.9312859918509488E-01};
    X[11] = vector<double> {-9.9429458548239924E-01, -9.7006049783542869E-01, -9.2695677218717398E-01, -8.6581257772030018E-01, -7.8781680597920811E-01, -6.9448726318668275E-01, -5.8764040350691160E-01, -4.6935583798675706E-01, -3.4193582089208424E-01, -2.0786042668822130E-01, -6.9739273319722211E-02,  6.9739273319722211E-02,  2.0786042668822130E-01,  3.4193582089208424E-01,  4.6935583798675706E-01,  5.8764040350691160E-01,  6.9448726318668275E-01,  7.8781680597920811E-01,  8.6581257772030018E-01,  9.2695677218717398E-01, 9.7006049783542869E-01, 9.9429458548239924E-01};
    X[12] = vector<double> {-9.9518721999702131E-01, -9.7472855597130947E-01, -9.3827455200273280E-01, -8.8641552700440096E-01, -8.2000198597390295E-01, -7.4012419157855436E-01, -6.4809365193697555E-01, -5.4542147138883956E-01, -4.3379350762604513E-01, -3.1504267969616340E-01, -1.9111886747361631E-01, -6.4056892862605630E-02,  6.4056892862605630E-02,  1.9111886747361631E-01,  3.1504267969616340E-01,  4.3379350762604513E-01,  5.4542147138883956E-01,  6.4809365193697555E-01,  7.4012419157855436E-01,  8.2000198597390295E-01, 8.8641552700440096E-01, 9.3827455200273280E-01, 9.7472855597130947E-01, 9.9518721999702131E-01};
    X[13] = vector<double> {-9.9588570114561692E-01, -9.7838544595647092E-01, -9.4715906666171423E-01, -9.0263786198430707E-01, -8.4544594278849805E-01, -7.7638594882067880E-01, -6.9642726041995728E-01, -6.0669229301761807E-01, -5.0844071482450570E-01, -4.0305175512348629E-01, -2.9200483948595690E-01, -1.7685882035689018E-01, -5.9230093429313208E-02,  5.9230093429313208E-02,  1.7685882035689018E-01,  2.9200483948595690E-01,  4.0305175512348629E-01,  5.0844071482450570E-01,  6.0669229301761807E-01,  6.9642726041995728E-01, 7.7638594882067880E-01, 8.4544594278849805E-01, 9.0263786198430707E-01, 9.4715906666171423E-01, 9.7838544595647092E-01, 9.9588570114561692E-01};
    X[14] = vector<double> {-9.9644249757395442E-01, -9.8130316537087281E-01, -9.5425928062893817E-01, -9.1563302639213207E-01, -8.6589252257439497E-01, -8.0564137091717913E-01, -7.3561087801363179E-01, -6.5665109403886501E-01, -5.6972047181140173E-01, -4.7587422495511827E-01, -3.7625151608907870E-01, -2.7206162763517810E-01, -1.6456928213338079E-01, -5.5079289884034266E-02,  5.5079289884034266E-02,  1.6456928213338079E-01,  2.7206162763517810E-01,  3.7625151608907870E-01,  4.7587422495511827E-01,  5.6972047181140173E-01, 6.5665109403886501E-01, 7.3561087801363179E-01, 8.0564137091717913E-01, 8.6589252257439497E-01, 9.1563302639213207E-01, 9.5425928062893817E-01, 9.8130316537087281E-01, 9.9644249757395442E-01};
    X[15] = vector<double> {-9.9689348407464951E-01, -9.8366812327974729E-01, -9.6002186496830755E-01, -9.2620004742927431E-01, -8.8256053579205263E-01, -8.2956576238276836E-01, -7.6777743210482619E-01, -6.9785049479331585E-01, -6.2052618298924289E-01, -5.3662414814201986E-01, -4.4703376953808915E-01, -3.5270472553087812E-01, -2.5463692616788985E-01, -1.5386991360858354E-01, -5.1471842555317698E-02,  5.1471842555317698E-02,  1.5386991360858354E-01,  2.5463692616788985E-01,  3.5270472553087812E-01,  4.4703376953808915E-01, 5.3662414814201986E-01, 6.2052618298924289E-01, 6.9785049479331585E-01, 7.6777743210482619E-01, 8.2956576238276836E-01, 8.8256053579205263E-01, 9.2620004742927431E-01, 9.6002186496830755E-01, 9.8366812327974729E-01, 9.9689348407464951E-01};
    X[16] = vector<double> {-9.9726386184948157E-01, -9.8561151154526838E-01, -9.6476225558750639E-01, -9.3490607593773967E-01, -8.9632115576605220E-01, -8.4936761373256997E-01, -7.9448379596794239E-01, -7.3218211874028971E-01, -6.6304426693021523E-01, -5.8771575724076230E-01, -5.0689990893222936E-01, -4.2135127613063533E-01, -3.3186860228212767E-01, -2.3928736225213706E-01, -1.4447196158279649E-01, -4.8307665687738310E-02,  4.8307665687738310E-02,  1.4447196158279649E-01,  2.3928736225213706E-01,  3.3186860228212767E-01, 4.2135127613063533E-01, 5.0689990893222936E-01, 5.8771575724076230E-01, 6.6304426693021523E-01, 7.3218211874028971E-01, 7.9448379596794239E-01, 8.4936761373256997E-01, 8.9632115576605220E-01, 9.3490607593773967E-01, 9.6476225558750639E-01, 9.8561151154526838E-01, 9.9726386184948157E-01};
    X[17] = vector<double> {-9.9757175379084195E-01, -9.8722781640630952E-01, -9.6870826253334430E-01, -9.4216239740510710E-01, -9.0780967771832455E-01, -8.6593463833456441E-01, -8.1688422790093362E-01, -7.6106487662987299E-01, -6.9893911321626290E-01, -6.3102172708052851E-01, -5.5787550066974667E-01, -4.8010654519032703E-01, -3.9835927775864594E-01, -3.1331108133946323E-01, -2.2566669161644948E-01, -1.3615235725918298E-01, -4.5509821953102533E-02,  4.5509821953102533E-02,  1.3615235725918298E-01,  2.2566669161644948E-01, 3.1331108133946323E-01, 3.9835927775864594E-01, 4.8010654519032703E-01, 5.5787550066974667E-01, 6.3102172708052851E-01, 6.9893911321626290E-01, 7.6106487662987299E-01, 8.1688422790093362E-01, 8.6593463833456441E-01, 9.0780967771832455E-01, 9.4216239740510710E-01, 9.6870826253334430E-01, 9.8722781640630952E-01, 9.9757175379084195E-01};
    X[18] = vector<double> {-9.9783046248408580E-01, -9.8858647890221230E-01, -9.7202769104969799E-01, -9.4827298439950758E-01, -9.1749777451565906E-01, -8.7992980089039707E-01, -8.3584716699247530E-01, -7.8557623013220657E-01, -7.2948917159355664E-01, -6.6800123658552102E-01, -6.0156765813598057E-01, -5.3068028592624517E-01, -4.5586394443342027E-01, -3.7767254711968923E-01, -2.9668499534402826E-01, -2.1350089231686559E-01, -1.2873610380938480E-01, -4.3018198473708608E-02,  4.3018198473708608E-02,  1.2873610380938480E-01, 2.1350089231686559E-01, 2.9668499534402826E-01, 3.7767254711968923E-01, 4.5586394443342027E-01, 5.3068028592624517E-01, 6.0156765813598057E-01, 6.6800123658552102E-01, 7.2948917159355664E-01, 7.8557623013220657E-01, 8.3584716699247530E-01, 8.7992980089039707E-01, 9.1749777451565906E-01, 9.4827298439950758E-01, 9.7202769104969799E-01, 9.8858647890221230E-01, 9.9783046248408580E-01};
    X[19] = vector<double> {-9.9804993053568758E-01, -9.8973945426638554E-01, -9.7484632859015352E-01, -9.5346633093352962E-01, -9.2574133204858433E-01, -8.9185573900463222E-01, -8.5203502193236214E-01, -8.0654416760531689E-01, -7.5568590375397071E-01, -6.9979868037918436E-01, -6.3925441582968168E-01, -5.7445602104780713E-01, -5.0583471792793111E-01, -4.3384716943237650E-01, -3.5897244047943500E-01, -2.8170880979016527E-01, -2.0257045389211670E-01, -1.2208402533786741E-01, -4.0785147904578239E-02,  4.0785147904578239E-02, 1.2208402533786741E-01, 2.0257045389211670E-01, 2.8170880979016527E-01, 3.5897244047943500E-01, 4.3384716943237650E-01, 5.0583471792793111E-01, 5.7445602104780713E-01, 6.3925441582968168E-01, 6.9979868037918436E-01, 7.5568590375397071E-01, 8.0654416760531689E-01, 8.5203502193236214E-01, 8.9185573900463222E-01, 9.2574133204858433E-01, 9.5346633093352962E-01, 9.7484632859015352E-01, 9.8973945426638554E-01, 9.9804993053568758E-01};
    X[20] = vector<double> {-9.9823770971055925E-01, -9.9072623869945708E-01, -9.7725994998377430E-01, -9.5791681921379168E-01, -9.3281280827867652E-01, -9.0209880696887434E-01, -8.6595950321225956E-01, -8.2461223083331170E-01, -7.7830565142651942E-01, -7.2731825518992710E-01, -6.7195668461417957E-01, -6.1255388966798030E-01, -5.4946712509512818E-01, -4.8307580168617870E-01, -4.1377920437160498E-01, -3.4199409082575849E-01, -2.6815218500725369E-01, -1.9269758070137111E-01, -1.1608407067525521E-01, -3.8772417506050816E-02, 3.8772417506050816E-02, 1.1608407067525521E-01, 1.9269758070137111E-01, 2.6815218500725369E-01, 3.4199409082575849E-01, 4.1377920437160498E-01, 4.8307580168617870E-01, 5.4946712509512818E-01, 6.1255388966798030E-01, 6.7195668461417957E-01, 7.2731825518992710E-01, 7.7830565142651942E-01, 8.2461223083331170E-01, 8.6595950321225956E-01, 9.0209880696887434E-01, 9.3281280827867652E-01, 9.5791681921379168E-01, 9.7725994998377430E-01, 9.9072623869945708E-01, 9.9823770971055925E-01};

    X[21] = vector<double> {-9.983996189900625e-01, -9.915772883408609e-01, -9.793425080637482e-01, -9.617593653382045e-01, -9.389235573549881e-01, -9.109597249041275e-01, -8.780205698121728e-01, -8.402859832618169e-01, -7.979620532554874e-01, -7.512799356894805e-01, -7.004945905561712e-01, -6.458833888692479e-01, -5.877445974851093e-01, -5.263957499311923e-01, -4.6217191207042196e-01, -3.9542385204297503e-01, -3.265161244654115e-01, -2.5582507934287907e-01, -1.8373680656485455e-01, -1.1064502720851986e-01, -3.694894316535177e-02, 3.694894316535177e-02, 1.1064502720851986e-01, 1.8373680656485455e-01, 2.5582507934287907e-01, 3.265161244654115e-01, 3.9542385204297503e-01, 4.6217191207042196e-01, 5.263957499311923e-01, 5.877445974851093e-01, 6.458833888692479e-01, 7.004945905561712e-01, 7.512799356894805e-01, 7.979620532554874e-01, 8.402859832618169e-01, 8.780205698121728e-01, 9.109597249041275e-01, 9.389235573549881e-01, 9.617593653382045e-01, 9.793425080637482e-01, 9.915772883408609e-01, 9.983996189900625e-01};
    X[22] = vector<double> {-9.985402006367742e-01, -9.923163921385159e-01, -9.81151833077914e-01, -9.650996504224931e-01, -9.442395091181941e-01, -9.186752599841758e-01, -8.885342382860432e-01, -8.539665950047104e-01, -8.15144539645135e-01, -7.722614792487559e-01, -7.25531053660717e-01, -6.751860706661224e-01, -6.214773459035758e-01, -5.646724531854708e-01, -5.050543913882023e-01, -4.429201745254115e-01, -3.7857935201470716e-01, -3.123524665027858e-01, -2.4456945692820126e-01, -1.7556801477551678e-01, -1.0569190170865325e-01, -3.5289236964135356e-02, 3.5289236964135356e-02, 1.0569190170865325e-01, 1.7556801477551678e-01, 2.4456945692820126e-01, 3.123524665027858e-01, 3.7857935201470716e-01, 4.429201745254115e-01, 5.050543913882023e-01, 5.646724531854708e-01, 6.214773459035758e-01, 6.751860706661224e-01, 7.25531053660717e-01, 7.722614792487559e-01, 8.15144539645135e-01, 8.539665950047104e-01, 8.885342382860432e-01, 9.186752599841758e-01, 9.442395091181941e-01, 9.650996504224931e-01, 9.81151833077914e-01, 9.923163921385159e-01, 9.985402006367742e-01};
    X[23] = vector<double> {-9.98663042133818e-01, -9.929623489061743e-01, -9.827336698041669e-01, -9.680213918539919e-01, -9.488923634460897e-01, -9.25433798806754e-01, -8.97752711533942e-01, -8.65975394866858e-01, -8.302468370660661e-01, -7.907300570752742e-01, -7.47605359615666e-01, -7.010695120204057e-01, -6.513348462019977e-01, -5.986282897127152e-01, -5.431903302618026e-01, -4.852739183881647e-01, -4.251433132828284e-01, -3.630728770209957e-01, -2.9934582270187005e-01, -2.3425292220626975e-01, -1.6809117946710353e-01, -1.0116247530558424e-01, -3.377219001605204e-02, 3.377219001605204e-02, 1.0116247530558424e-01, 1.6809117946710353e-01, 2.3425292220626975e-01, 2.9934582270187005e-01, 3.630728770209957e-01, 4.251433132828284e-01, 4.852739183881647e-01, 5.431903302618026e-01, 5.986282897127152e-01, 6.513348462019977e-01, 7.010695120204057e-01, 7.47605359615666e-01, 7.907300570752742e-01, 8.302468370660661e-01, 8.65975394866858e-01, 8.97752711533942e-01, 9.25433798806754e-01, 9.488923634460897e-01, 9.680213918539919e-01, 9.827336698041669e-01, 9.929623489061743e-01, 9.98663042133818e-01};
    X[24] = vector<double> {-9.987710072524261e-01, -9.935301722663508e-01, -9.841245837228269e-01, -9.705915925462473e-01, -9.529877031604309e-01, -9.313866907065543e-01, -9.058791367155696e-01, -8.765720202742479e-01, -8.435882616243935e-01, -8.070662040294426e-01, -7.671590325157404e-01, -7.240341309238146e-01, -6.778723796326639e-01, -6.288673967765136e-01, -5.772247260839727e-01, -5.23160974722233e-01, -4.669029047509584e-01, -4.086864819907167e-01, -3.4875588629216075e-01, -2.8736248735545555e-01, -2.2476379039468905e-01, -1.612223560688917e-01, -9.70046992094627e-02, -3.238017096286937e-02, 3.238017096286937e-02, 9.70046992094627e-02, 1.612223560688917e-01, 2.2476379039468905e-01, 2.8736248735545555e-01, 3.4875588629216075e-01, 4.086864819907167e-01, 4.669029047509584e-01, 5.23160974722233e-01, 5.772247260839727e-01, 6.288673967765136e-01, 6.778723796326639e-01, 7.240341309238146e-01, 7.671590325157404e-01, 8.070662040294426e-01, 8.435882616243935e-01, 8.765720202742479e-01, 9.058791367155696e-01, 9.313866907065543e-01, 9.529877031604309e-01, 9.705915925462473e-01, 9.841245837228269e-01, 9.935301722663508e-01, 9.987710072524261e-01};
    X[25] = vector<double> {-9.98866404420071e-01, -9.940319694320907e-01, -9.853540840480058e-01, -9.72864385106692e-01, -9.566109552428079e-01, -9.36656618944878e-01, -9.130785566557919e-01, -8.859679795236131e-01, -8.554297694299461e-01, -8.21582070859336e-01, -7.845558329003992e-01, -7.444943022260686e-01, -7.015524687068222e-01, -6.558964656854394e-01, -6.077029271849502e-01, -5.5715830451465e-01, -5.044581449074642e-01, -4.4980633497403877e-01, -3.9341431189756515e-01, -3.3550024541943735e-01, -2.7628819377953195e-01, -2.1600723687604176e-01, -1.548905899981459e-01, -9.317470156008613e-02, -3.1098338327188876e-02, 3.1098338327188876e-02, 9.317470156008613e-02, 1.548905899981459e-01, 2.1600723687604176e-01, 2.7628819377953195e-01, 3.3550024541943735e-01, 3.9341431189756515e-01, 4.4980633497403877e-01, 5.044581449074642e-01, 5.5715830451465e-01, 6.077029271849502e-01, 6.558964656854394e-01, 7.015524687068222e-01, 7.444943022260686e-01, 7.845558329003992e-01, 8.21582070859336e-01, 8.554297694299461e-01, 8.859679795236131e-01, 9.130785566557919e-01, 9.36656618944878e-01, 9.566109552428079e-01, 9.72864385106692e-01, 9.853540840480058e-01, 9.940319694320907e-01, 9.98866404420071e-01};


    vector<double> cumsum_w(nGG, W[nGG][nGG]);
    for(int i=1; i<nGG; i++)
        cumsum_w[i] = cumsum_w[i-1] + W[nGG][i+nGG];

    Ft_pts.resize(nGGa);                            // \tilde{F} grid
    Ft_pts[0] = Fmin;
    for(int i=1; i<nGGa; i++)
        Ft_pts[i] = Fmin + (Fmax-Fmin)*cumsum_w[i-1];   

    F_pts.resize(nGG);
    for(auto i=nGG; i < X[nGG].size(); i++)
        F_pts[i-nGG] = Fmin + X[nGG][i]*(Fmax - Fmin);    // F grid (vals bet. \tild{F} pnts)

}

////////////////////////////////////////////////////////////////
///
/// Read the albdf table for CO2, CO, and H2O.
/// Separate files are given for various pressures.
/// Interpolate the files to the desired pressure: P.
/// Reshape from the given 1D to the 3D (CO/CO2) or 4D (H2O) tables.
///
////////////////////////////////////////////////////////////////

void rad_rcslw::set_Falbdf_CO2_CO_H2O_at_P(){

    double P1;
    double P2;
    double f;

    if(P < P_table[0] || P > P_table.back() ) {
        cout << "Pressure = " << P << " atm is out of range\n";
        exit(0);
    }
    else if(P==P_table[0]){
        P1 = P_table[0];
        P2 = P_table[1];
    }
    else if(P==P_table.back()){
        P1 = P_table[nP-2];
        P2 = P_table[nP-1];
    }
    else{
        int i;
        for(i=0; P_table[i] < P; i++)
            P1 = P_table[i];
        P2 = P_table[i];
    }

    f  = (P-P1)/(P2-P1);               // Interpolation factor

    string Pres_1 = to_string(P1);
    string Pres_2 = to_string(P2);
    Pres_1.erase(Pres_1.find_last_not_of('0')+1, string::npos);
    Pres_2.erase(Pres_2.find_last_not_of('0')+2, string::npos);
    if(Pres_1.back() == '.') Pres_1.push_back('0');
    if(Pres_2.back() == '.') Pres_2.push_back('0');
    replace(Pres_1.begin(), Pres_1.end(), '.', '_');
    replace(Pres_2.begin(), Pres_2.end(), '.', '_');

    //--------------------- CO2
    
    string CO2_file1 = string(GETSTRINGEDVALUE(RCSLW_DATA_DIR)) + "/co2_p" + Pres_1 + ".bin";
    string CO2_file2 = string(GETSTRINGEDVALUE(RCSLW_DATA_DIR)) + "/co2_p" + Pres_2 + ".bin";
    //string CO2_file1 = "ALBDF_Tables/co2_p" + Pres_1 + ".bin";
    //string CO2_file2 = "ALBDF_Tables/co2_p" + Pres_2 + ".bin";

    vector<double> CO2_F1(nTg*nTb*nC);
    vector<double> CO2_F2(nTg*nTb*nC);
    Falbdf_CO2.resize(nTg*nTb*nC);

    get_FI_albdf_tables(CO2_file1, nTg, nTb, nC, CO2_F1);
    get_FI_albdf_tables(CO2_file2, nTg, nTb, nC, CO2_F2);

    for(auto ind=0; ind < Falbdf_CO2.size(); ++ind)
        Falbdf_CO2[ind] = CO2_F1[ind]*(1-f) + CO2_F2[ind]*f;

    //---------------------CO

    string CO_file1 = string(GETSTRINGEDVALUE(RCSLW_DATA_DIR)) + "/co_p" + Pres_1 + ".bin";
    string CO_file2 = string(GETSTRINGEDVALUE(RCSLW_DATA_DIR)) + "/co_p" + Pres_2 + ".bin";
    //string CO_file1 = "ALBDF_Tables/co_p" + Pres_1 + ".bin";
    //string CO_file2 = "ALBDF_Tables/co_p" + Pres_2 + ".bin";

    vector<double> CO_F1(nTg*nTb*nC);
    vector<double> CO_F2(nTg*nTb*nC);
    Falbdf_CO.resize(nTg*nTb*nC);

    get_FI_albdf_tables(CO_file1, nTg, nTb, nC, CO_F1);
    get_FI_albdf_tables(CO_file2, nTg, nTb, nC, CO_F2);

    for(auto ind=0; ind < Falbdf_CO.size(); ++ind)
        Falbdf_CO[ind] = CO_F1[ind]*(1-f) + CO_F2[ind]*f;

    //---------------------H2O
   
    string H2O_file1 = string(GETSTRINGEDVALUE(RCSLW_DATA_DIR)) + "/h2o_p" + Pres_1 + ".bin";
    string H2O_file2 = string(GETSTRINGEDVALUE(RCSLW_DATA_DIR)) + "/h2o_p" + Pres_2 + ".bin";
    //string H2O_file1 = "ALBDF_Tables/h2o_p" + Pres_1 + ".bin";
    //string H2O_file2 = "ALBDF_Tables/h2o_p" + Pres_2 + ".bin";

    vector<double> H2O_F1(ny_H2O*nTg*nTb*nC); 
    vector<double> H2O_F2(ny_H2O*nTg*nTb*nC); 
    Falbdf_H2O.resize(ny_H2O*nTg*nTb*nC);

    get_FI_albdf_tables(H2O_file1, ny_H2O, nTg, nTb, nC, H2O_F1);
    get_FI_albdf_tables(H2O_file2, ny_H2O, nTg, nTb, nC, H2O_F2);

    for(auto ind=0; ind < Falbdf_H2O.size(); ++ind)
        Falbdf_H2O[ind] = H2O_F1[ind]*(1-f) + H2O_F2[ind]*f;
  
}
    
/////////////////////////////////////////////////////////////////////////
///
/// Import ALBDF files to array; 3D files
/// @param Ptable_file_name \input   file for given species for given pressure
/// @param nx               \input   number of points in first dimension
/// @param ny               \input   number of points in second dimension
/// @param nz               \input   number of points in third dimension
/// @param myarray          \output  array read from file
///
/////////////////////////////////////////////////////////////////////////

void rad_rcslw::get_FI_albdf_tables(const string Ptable_file_name,
                                    const int nx, const int ny, const int nz,
                                    vector<double> &myarray){

    ifstream ifile;
    ifile.open(Ptable_file_name, ios::in | ios::binary);
    if(!ifile){
        cout << endl << "error opening file: " << Ptable_file_name << endl;
        exit(0);
    }
    float data;
    int ntot = nx*ny*nz;
    for(int ind=0; ind<ntot; ++ind){
        ifile.read((char *) &data, sizeof(data));
        myarray[ind] = static_cast<double>(data);
    }
    ifile.close();
    
}

/////////////////////////////////////////////////////////////////////////
///
/// Import ALBDF files to array; 4D files
/// @param Ptable_file_name \input   file for given species for given pressure
/// @param nx               \input   number of points in first dimension
/// @param ny               \input   number of points in second dimension
/// @param nz               \input   number of points in third dimension
/// @param nw               \input   number of points in fourth dimension
/// @param myarray          \output  array read from file
///
/////////////////////////////////////////////////////////////////////////

void rad_rcslw::get_FI_albdf_tables(const string Ptable_file_name,
                                    const int nx, const int ny, const int nz, const int nw,
                                    vector<double> &myarray){
 
    ifstream ifile;
    ifile.open(Ptable_file_name, ios::in | ios::binary);
    if(!ifile){
        cout << endl << "error opening file: " << Ptable_file_name << endl;
        exit(0);
    }
    float data;
    int ntot = nx*ny*nz*nw;
    for(int ind=0; ind<ntot; ++ind){
        ifile.read((char *) &data, sizeof(data));
        myarray[ind] = static_cast<double>(data);
    }
    ifile.close();
    
}
////////////////////////////////////////////////////////////////
/**
 Get soot albdf F
 @param C       \input    cross section (m2/mol)
 @param Tg      \input    gas temperature
 @param Tb      \input    black temperature
 @param fvsoot  \input    soot volume fraction
 @return the albdf function F for soot
 Note: This comes from Solovjov and Webb, J. Heat Transfer, <a href="https://doi.org/10.1115/1.1350824" target="_blank">123(3):450-457</a> (2001),
       and from Chang and Rhee (eq. 5), Int. Comm. Heat and Mass Transfer, <a href="https://doi.org/10.1016/0735-1933(84)90051-4" target="_blank">11(5):451-455</a> (1984).

 For computing csoot, see Radiative Heat Transfer, 3rd edition, by Modest, pages 424-425:
 \verbatim
     ksoot = 1.23, 1.95,  1.30,    0.92,  0.71, for
     nsoot = 2.21, 2.63,  2.19,    1.89,  2.31,
             Lee,  Stull, Dalzell, Chang, Felske, respectively
 \endverbatim
  *  See also Williams, Shaddix et al. Int. J. Heat and Mass Transfer <a href="https://www.sciencedirect.com/science/article/pii/S0017931006004893" target="_blank">50:1616-1630</a> (2007), 
 \verbatim
     ksoot = 1.03,    0.56,                0.43,         0.89,     0.80 for
     nsoot = 1.75,    1.57,                1.90,         1.99,     1.55
             Shaddix, Dalzell and Sarofim, Lee and Tien, Krishnan, Mountain and Mulholland, respectively
 \endverbatim
  *  Note, these give Planck Mean absorption coefficients of \f$(3.72*c_{soot}/C_2)*fv*T\f$, where \f$C_2 = 0.014388\, m\cdot K\f$.
 \verbatim
     Hence, (3.72*csoot/C2) = 1361, 1141, 1423, 1476, 835.0 for Lee, Stull, ...
                            = 1817, 1265, 744.2, 1319, 1785 for Shaddix, Dalzell and Sarofim, ...
 \endverbatim
  * Shaddix's constants give \f$c_{soot} = 7.03\f$, which is the same as the value of 7 presented in Solovjov 2001 for 
  *   Eq. 16 proposed by Hottel and Sarofim in their 1967 textbook Radiative Transfer.
**/
////////////////////////////////////////////////////////////////

double rad_rcslw::F_albdf_soot(const double C, const double Tg, const double Tb, const double fvsoot){

    if (fvsoot < 1E-12)
        return 1.0;
    
    double ksoot = 1.03;                // real part of complex refractive index
    double nsoot = 1.75;                // imag part of complex refractive index
    double csoot = 36*M_PI*nsoot*ksoot/(pow((nsoot*nsoot - ksoot*ksoot + 2),2) + 4*pow(nsoot*ksoot,2));

    double hCokb = 0.01438777354;      // h*Co/kb = C2 = PlancK*lightSpeed/Boltzmann = 6.62607004E-34 m2*kg/s * 299792458 m/s / 1.38064852E-23 m2*kg/s2*K
    
    double Nconc = P*101325/8.31446/Tg;      // mol/m3

    double x = hCokb*C*Nconc / (csoot*fvsoot*Tb);
    double sum = 0.0;
    double nx;
    for(int n=1; n<=3; n++){          // sum from n=1 to oo, but n=1 to 3 is enough
        nx = n*x;
        sum += exp(-nx)/pow(n,4)*(nx*(nx*(nx+3)+6)+6);
    }
    return 1.0 - 15.0/pow(M_PI,4)*sum;
}

